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We discuss the linear and two-photon spectroscopic selection rules for spin-singlet excitons in monolayer 
transition-metal dichalcogenides. Our microscopic formalism combines a fully ^-dependent few-orbital band 
structure with a many-body Bethe-Salpeter equation treatment of the electron-hole interaction, using a model 
dielectric function. We show analytically and numerically that the single-particle, valley-dependent selection 
rules are preserved in the presence of excitonic effects. Furthermore, we definitively demonstrate that the bright 
(one-photon allowed) excitons have 5-type azimuthal symmetry and that dark p -type excitons can be probed via 
two-photon spectroscopy. The screened Coulomb interaction in these materials substantially deviates from the 
form; this breaks the “accidental” angular momentum degeneracy in the exciton spectrum, such that the 
2 p exciton has a lower energy than the 2s exciton by at least 50 meV. We compare our calculated two-photon 
absorption spectra to recent experimental measurements. 


I. INTRODUCTION 

The transition-metal dichalcogenides (TMDCs) are a fam¬ 
ily of layered semiconducting crystals that includes M 0 S 2 , 
MoSe 2 , WS 2 , and WSe 2 . Isolated monolayers of TMDCs 
have been recently investigated for two major reasons. First, 
the emergent direct band-gap occurs at the corners of the 
hexagonal Brillouin zone (so-called Valleys’) [1,2] and the 
nearby band structure topology leads to valley-dependent op¬ 
tical selection rules [3-5]. Second, the carrier confinement 
and reduced dielectric screening leads to large many-body ef¬ 
fects, such as the formation of strongly bound excitons [6-10], 
trions [6, 11-13], and biexcitons [14] with very large binding 
energies. A unified understanding of the optical properties 
must treat both of these aspects on equal footing, and signifi¬ 
cant effort is now being focused on investigating the detailed 
spectroscopy of excitons in monolayer TMDCs. 

In the ongoing effort to understand excitons in these mate¬ 
rials, multiple spectroscopic techniques have been employed, 
including reflectance [9, 10, 15], photoluminescence excita¬ 
tion spectroscopy [16] scanning tunneling spectroscopy [17, 
18], and two-photon luminescence [9, 19, 20]. A rigorous 
knowledge of the spectroscopic selection rules for excitons 
in monolayer TMDCs is crucial for the proper interpretation 
of these and future experiments. In this paper, we develop a 
model-based framework which is sufficiently detailed to pro¬ 
vide quantitative results, but also sufficiently simple to al¬ 
low precise statements about symmetry-determined selection 
rules. We describe the connection to our previous work based 
on an effective mass theory of excitons [6], and identify the 
key microscopic physical factors that determine the properties 
of excitons and their interaction with photons. We also pro¬ 
vide the first theoretical treatment of two-photon absorption 
in monolayer TMDCs. 

The outline of the paper is as follows. In Sec. II we will 
discuss simple microscopic models of the single-particle band 
structure in monolayer TMDCs, and in particular we will an¬ 
alyze the transition matrix elements which completely deter¬ 
mine the independent-electron absorption and partially deter¬ 


mine the excitonic absorption. We will then in Sec. Ill ana¬ 
lyze the linear optical properties and present selection rules, 
both in the absence and presence of exciton effects, defini¬ 
tively finding that 5-type excitons are optically bright. Lastly, 
in Sec. IV we will calculate the two-photon absorption sig¬ 
nal which will be shown to probe p-type excitons and we will 
discuss some of the implications for recent experiments. We 
conclude in Sec. V, and make connection to other recent theo¬ 
retical works. We note that a preliminary version of this work 
appeared in Ref. [21]. 


II. SINGLE-PARTICLE BAND STRUCTURE 


We will consider two models for the single-particle band 
structure. First, we will consider a widely used long- 
wavelength, two-band model [3]. In particular, this minimal 
model allows for a largely analytical treatment, which exposes 
many of the subtleties of the theory, including selection rules 
and exciton effects. Second, we will use a recently presented 
nonlinear three-band model [22], which requires a numerical 
treatment but captures higher-order effects. This also ensures 
that our conclusions are generally valid and not specifically 
dependent on the simplified two-band picture. For simplicity 
we will henceforth neglect spin-orbit coupling, though it can 
be straightforwardly included in the single-particle descrip¬ 
tion [3, 22, 23]. Specifically, in all models of the band struc¬ 
ture, the spin projection 5 Z is still a good quantum number in 
the presence of spin-orbit coupling. In this sense, the follow¬ 
ing discussion applies to the A-exciton (and not the 5-exciton) 
and conventional factors of two for spin will not appear. At 
this level of theory, the formalism for the iLexciton is identi¬ 
cal, and its contribution is simply shifted to higher energies. 
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A. Two-band model 


The first model considered has the form of a conventional 
two-band, massive Dirac Hamiltonian, 


Hr(k) = 


E g 12 at(jq x - iq y ) \ 
at(rq x + iq y ) -E g / 2 ) 


( 1 ) 


The variable r = ±1 indexes the two “valleys,” known as the K 
and K' (or K and -K) points, which occur at alternating cor¬ 
ners of the hexagonal first Brillouin zone. The Hamiltonian 
has been linearized in the wavevector difference with respect 
to the nearest K point, i.e. q - k - K. This is a gapped 
version of the conventional graphene Hamiltonian [24]. In 
graphene, the spinor basis corresponds to carbon p z orbitals 
on the two distinct sublattices; in the TMDCs, the basis corre¬ 
sponds to the transition-metal \d z i) = \<p c ) orbital and the metal 
symmetry-adapted \d x 2 _ y i) + ir\d xy ) = \(j) T v ) orbital. The above 
Hamiltonian was first used for TMDCs by Xiao et al. [3] who 
predicted optical selection rules leading to spin-valley cou¬ 
pling. Such spin-valley coupling was quickly confirmed ex¬ 
perimentally, by monitoring the circular polarization of pho¬ 
toluminescence [4, 5]. 

The eigenvalues of the two-band Hamiltonian are 

Ec,v(k) = ± I y[Ej + Mfltqf = ±£(k) (2) 

and the eigenvectors are 

K, k ) = -T [ VMfc) \<Pc) + VMfcy' T ^l^>] (3a) 

\K,k > = -j= [- A /al(k)\(p c ) + \<t>l)\ ■ (3b) 

where a+(k) = 1 ± E g /[2s(k)] and tan0^ = q y /q x . The rel¬ 
ative phase appearing within each eigenvector is associated 
with an electronic “chirality” (related to Berry’s phase), which 
is well-known in graphene [24-26]. Note that the overall 
phase of each eigenvector is arbitrary, and the phase conven¬ 
tion chosen here is such that the first element of each eigen¬ 
vector is purely real. 


B. Three-band model 

A more detailed Hamiltonian - using three bands derived 
from the transition-metal \d z i), \d xy ), and \d x 2 _ y i) atomic or¬ 
bitals - was given recently by Liu et al [22]. The form of 
the matrix elements and material-specific parameters can be 
found in Ref. [22]. We note than in addition to using three 
bands instead of two, this Hamiltonian has not been linearized 
with respect to wavevector near the K and K' points, which 
gives a more accurate description throughout the entire Bril¬ 
louin zone. While it cannot be so easily diagonalized analyt¬ 
ically, the Hamiltonian can be straightforwardly diagonalized 
numerically. For phase consistency in later calculations, we 
enforce the same phase convention as for the two-band eigen¬ 
vectors, i.e. that the first element of each eigenvector is purely 
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FIG. 1. Single-particle band structure of MoS 2 predicted by a lin¬ 
earized two-band model (blue solid) and a non-linear three-band 
model (red dashed) compared to first-principles density functional 
theory with the local density approximation (DFT, solid black). 


real, which is sufficient to ensure continuity in k- space. In 
Fig. 1, the band structure predicted by these two models is 
compared to the band structure calculated by density func¬ 
tional theory with the local density approximation. 


C. Transition matrix elements 

An analysis of optical selection rules requires the momen¬ 
tum matrix elements between single-particle states. In the 
present model Hamiltonians, the momentum matrix elements 
normal to the layer are zero by symmetry. Here we focus on 
the momentum in the plane. By using the commutation rela¬ 
tion p = ( -im/ti)[r , //], we can write these momentum matrix 
elements as 


P vc (k)= —(^, k \[r,H]\^ k ) 

n (4) 

m v y 

= ~7~ (E C ,k ~ Ev,k) (0v,fc|V k\lfSc,k) 
n 

where we have used the k-space representation of the position 
operator, r = /V*,. We can now use a generalized Feynman- 
Hellman theorem to write this as 

YYl 

P vc (k ) = -<^v,fc|V fc ^(fc)|^ ffc > (5) 

n 

(note that this expression neglects the on-site, intra-atomic 
contribution [27], however this vanishes here for d-d transi¬ 
tions). For the simple two-band Hamiltonian, this gives 

T7 _ / 0 - iy) 

VkH(k) -{at(rx + iy ) 0 

The appropriate matrix element can then be taken between the 
conduction and valence band eigenstates of the Hamiltonian, 
yielding a transition dipole vector P vc (k ) with linear v- and 
y-polarization components 

P v x (k) = T ~^~~ ^r+(k)e~ lT( ^ k - a-(k)e lT( ^ k J, (7) 

P v y c ( k ) = i^j- [a + (k)e- iT *« + a-{k)e]. (8) 
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FIG. 2. Valence to lowest conduction band momentum matrix elements for a linearized two-band model (top) and a non-linear three-band 
model (bottom). Blue is positive, red is negative, and white is zero. The results are qualitatively very similar in the immediate vicinity of the 
K and K' points, but differ elsewhere in the Brillouin zone. 


The same procedure can be done for the three-band Hamilto¬ 
nian, by taking the gradient and calculating (numerically) the 
appropriate matrix element between conduction and valence 
bands. A comparison of the real and imaginary parts of the 
x- and y-components of the two different models of the band 
structure is shown in Fig. 2 throughout the entire first Brillouin 
zone. 

Valley-dependent selection rules have been shown to arise 
specifically for the case of circularly polarized light [3]. For 
circular polarizations, the above expressions can be combined 
to give, in the two-band case, 


PT(k) = -h [p v x c (k) ± ip; c (k )] 


V 2 

mat 


1 + T- 




V2 fi \ 2 e(k) 

leading to the valley-dependent intensities, 


(9) 


shown in Fig. 3. Note that while the matrix elements them¬ 
selves have an ambiguity in the phase (i.e. they are not ob¬ 
servable), the squared matrix elements are completely inde¬ 
pendent of the phase convention. In Sec. IIIB, we will show 
how the nodal structure (p-type symmetry) of the momentum 
matrix elements is canceled, leading to bright s-type excitons 
which still respect the valley selectivity. 


III. LINEAR OPTICAL PROPERTIES AND SELECTION 
RULES 


In general, the transition probability per unit time is given 
by 

W(to) = yTj ]Vif1 2 6{Ef (H) 


\Pl c (k)\ 2 


m 2 a 2 t 2 
2 h 2 


1 + T 


2 e(k) 


( 10 ) 


Near the K and K' points, 2 s(k) —> E g , such that P v ±{k) oc 

(1 + r)e ±l(f)k and |P^ c (/c)| 2 oc (1 + r) 2 , i.e. circular polarization 
can selectively excite electrons at the K or K' point. For ex¬ 
ample, right-handed circular polarization, P v f(k), selectively 
excites at the K (r = +1) point. Again, this analysis can be 
carried out numerically for the three-band model. A compari¬ 
son of the the selection rules, |P+ c (/c)| , for the two models is 


where Vif is the matrix element which couples the initial and 
final states with energies £/ and E F . For the linear (one- 
photon) absorption, we have V = ( eA/mc)X • p, where A is the 
vector potential and A is the polarization. Within the presently 
considered model Hamiltonians, symmetry excludes coupling 
to photons with electric vector polarized perpendicular to the 
plane of the monolayer. Here we explicitly consider the case 
with electric vector polarized in the plane. We will evalu¬ 
ate this expression first in the independent particle picture and 
then in the presence of excitonic effects. 
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FIG. 3. Valence to lowest conduction band momentum matrix el¬ 
ements squared, for circular polarization, for a linearized two-band 
model (top) and a non-linear three-band model (bottom). Black is 
positive and white is zero. The results are qualitatively very similar 
in the immediate vicinity of the K and K' points, but differ elsewhere 
in the Brillouin zone. 


A. Independent particle absorption 

For an uncorrelated initial ground state |0) and an uncorre¬ 
lated final excited state k c v ^ |0), it is simple to show 


we can carry out the integration in Eq. (15) semi-analytically. 
Considering only one valley (say r = +1), we can change to 
polar coordinates about the K point, 

r\ 2 r OO 

4(w) = —7—7 k\P v _ c (k)\ 2 6(2e(k) - Tuo)dk. (16) 

m z co z Jo 

Note that by integrating out to infinity, we are incurring an 
error at large wavevectors (energies). Since the dispersion re¬ 
lation is monotonic, we can change variables, kdk = sets/a 2 t 2 , 
and use the squared matrix element from above to find 




f 


ds6(2s - Eo)£\ 1 + 


2s 


Tie 2 I E g 

= i + 


S(2s - fico) 


(17) 


Accounting for the other valley, £2 (oj) = £:f(tu) + 
yields 


Tie 


si(cS) = -rr-Oifuo - E g ) 
ATKjJ 





f 

1 + 




(18) 


At energies just above the gap, the dielectric function is like 
that of a conventional 2D semiconductor, i.e. oj 2 £2(oj) = 
const, but at higher energies it behaves like graphene (due to 
the linear dispersion), i.e. oj£ 2 (oj) = const. However, the lin¬ 
ear dispersion is unrealistic for TMDCs, as can be seen in the 
full band structure (Fig. 1). 


Vi F = —<0|A • pci fc c Vj fc|0) = — A P vc (k), (12) 
me ’ me 

B. Exciton absorption and the Bethe-Salpeter equation 


E f -Ej = EJk) - Ey(k), 


(13) 


and therefore 

2n 1 eA > 2 


n \mc) 


cv,k ( 14 ) 

x S(E c (k) - E v (k) - hoj). 


The imaginary part of the dielectric function [28] follows 
as [29] 


si (oj) = 


\k e 


2„2 


?x 


d z k 


\X • P vc {k)\ 


BZ ( 2 Tl) 2 

x S(E c (k) - E v (k) - hu). 


(15) 


where we have taken the infinite-system limit. Let us specif¬ 
ically consider the linearized two-band model with right- 
handed circular polarization, A • P vc (k) = P v _2(k), for which 


We now consider the spin-singlet optical properties includ¬ 
ing the excitonic effects arising from the strong electron-hole 
interaction. The correlated excited states within the single¬ 
excitation approximation can be written as 

BZ 

w =e E A - (fe) A c ^ |0 >’ < 19 ) 

k vc 

where |0) is again an uncorrelated (determinental) ground 
state. This form for the excited state wavefunction underlies 
the time-dependent Hartree-Fock and Bethe-Salpeter equation 
(BSE) formalisms; here we will pursue the latter, which is a 
many-body perturbative theory in the screened two-particle 
interaction. For a periodic crystal exciton wavefunction, 
Eq. (19), the BSE is an eigenvalue problem [30] for the ex¬ 
citon energy E x , 
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E x A* c (k) = (E Ctk - E v , 


1 


BZ 

k)A* c (k) + — E 

k' v',c' 


(^v,k^c,k\K eh \i// V ',h'^c',k')A* c ,(k'). 


( 20 ) 


The electron-hole interaction kernel K eh is the sum of a frequency-dependent screened Coulomb interaction and an unscreened 
exchange interaction [30], 


(Mc,k\K eh ’ d \^,k^ c , M ) = - Jd d r J 

('/' v ,k>J'c,k\K eh ’ x \<J/ V ',k’t/'c’,k’) = J' d d r J" d d r'i//l k (r)i// Vtk (r)\r-r'\~ l i// c ^ k -(r')iJ/* v , k ,(r'). 


(21a) 

(21b) 


If we (i) neglect the exchange interaction, (ii) neglect the frequency-dependence and local-field effects of the screened direct 
interaction, i.e. W(r , r', cS) ~ W(r - r', = 0), and (iii) make a ‘zero differential overlap’ approximation for the atomic orbitals, 

we find 


(ftv,klfrc,k\K ^c',k ') ~ {^c,kVkc',k f )(l/V,fc r \tffv,k)W(k k ). 


( 22 ) 


In the above, we have neglected the possible orbital structure to the screened interaction Wij(r - r'). 


At this point, we wish to emphasize that the orbital overlap 
prefactor in the screened interaction is crucially important. As 
an explicit example, in the two-band picture, we have 


WlkK*) = \[y/a + (k)a + (k') 

+ VQ'-(fe)Q'-(fe , )e"' T(0fc " 0fe ' ) ], 

Kk'Kk) = \\\/«-(k’)aJk) 


(23a) 


(23b) 


Coulomb interaction that weakly break certain degeneracies 
(see below). In this case, the spectrum of Eqs. (25) and (26) is 
no longer identical to that of the BSE with the screened inter¬ 
action given in Eq. (22). 

It remains to be shown whether the exciton wavefunctions 
of the original problem, as described by the BSE (25), have 
the same selection rules or the same spatial symmetries as the 
wavefunction of the real-space Wannier equation (26). To an¬ 
alyze the spatial symmetries, we can calculate the real-space 
wavefunction corresponding to the solution of the BSE, with 
the hole position fixed at the origin. We find 


As before, near the K and K' points, 2 s(k) —> E g , [i.e. 
a+(k ) ~ 1 and or_(fe) « 0], and in this limit, 

WckKk') * 1 (24a) 

<Kk'Kk>^ e ^ k ~ M - ( 24b > 

The BSE, Eq. (20), then yields a Wannier-//^, two-band pic¬ 
ture with an unusual phase factor in the screened interaction, 

Ex Alik) = (E c , k - E v , k )Al c (k) 

1 BZ 1251 

--Yw(k- k'y T ^ k ~ M A^ c {k'). 

Multiplying through by e~ lT(f>k gives a conventional Wannier 
equation for the pseudo-wavefunction A*(fe) = e~ lT(f)h A z c (k). 
If the bands can be approximated as parabolic, this means that 
the energy spectrum of the BSE is identical to that of a corre¬ 
sponding real-space Wannier equation with a screened inter¬ 
action W(r ), as we have employed in previous work [6, 10], 


Vx(r e ,r h = 0) = 2^(fc¥ c , fc (r e )^ ih ( 0) 

* (27) 

* }]Al c (k)e- iT ^e ik r ‘ =Al c (r e ), 

k 

demonstrating that the wavefunction which solves the real- 
space Eq. (26) is indeed (approximately) the same as the real- 
space BSE wavefunction. At a less approximate level, the 
spatial symmetries (s, p , d , etc.) will be identical. This is one 
of the main conclusions of this work. 

To determine the selection rules, we now consider the op¬ 
tical absorption in the presence of correlated excitonic ef¬ 
fects. Assuming as before an uncorrelated initial (ground) 
state |/) = |0), but now using a Wannier-like final exciton state 
\X) as in Eq. (19) gives 

</|A • p|X> = J A*ik)\ ■ P vc (k), (28) 

k 

which leads to the dielectric function 


47T 2 e 2 y 
m 2 oj 2 2j 

However, as explained in a recent work by Srivastava and 
Imamoglu [31], systematically continuing the expansion of Recall that for right-handed circular polarization, the momen- 

Eqs. (24) for small k - k' leads to additional terms in the turn matrix element near the K' (r = -1) point is nearly zero 


2_j K(k)A ■ P vc (k) 

k 


8(h(jj-E x ). (29) 


[E x -E g ]A x vc {r).= 


-J_V 2 
2H r 


W(r) 


Kir ) (26) 


eiiio) 
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FIG. 4. Imaginary part of the dielectric function for MoS 2 calculated 
in the presence of excitonic effects. The band gap has been rigidly 
increased to 2.41 eV such that the Is exciton peak occurs near 2.0 
eV (spin-orbit splitting into A and B peaks is neglected, as described 
in the text). A Gaussian broadening of 50 meV (FWHM) has been 
applied to all peaks. 


and near the K (r = 1) point it is given by 

A • P vc (k) = P v _ c (k ) « EE P 0 e~^ h . (30) 

n 

Therefore we can restrict our attention to k near K , which 
gives 

</|A • p\X) = P 0 J] Al(k)e-^ = P 0 Al(r = 0) (31) 

k~K 


and therefore 


£2 (oj) = 


47r 2 e 2 Po 

m 2 (o 2 


Y J \A vc (r = 0)\ 2 6(huj-E x ), 

X 


(32) 


which is just the usual Elliott formula for the excitonic absorp¬ 
tion [32]. In particular, the selection rules are conventional in 
that they are determined by the behavior of the wavefunction 
at the origin in real-space, leading to bright states with s-type 
azimuthal symmetry. We emphasize that the phase factor ap¬ 
pearing in the momentum matrix element is essentially can¬ 
celled by the conjugate phase factor in the exciton envelope 
wavefunction, which itself originates from the change of ba¬ 
sis in the screened interaction, Eq. (22). Therefore, not only 
can the excitons be labeled in analogy with the hydrogen se¬ 
ries in terms of their spatial symmetries but, to lowest order, 
they also obey identical selection rules. This is the second 
main conclusion of this work. 

As usual, the same analysis cannot be done analytically 
on the three-band model, but it can be straightforwardly car¬ 
ried out numerically. The dielectric function calculated via 
Eq. (29) for the two considered band structure models of 
M 0 S 2 is plotted in Fig. 4; in particular, the orbital overlaps 
in Eq. (22) are calculated numerically, without the approxi¬ 
mation given in Eqs. (24). As described in Refs. [6, 10], the 
screened interaction used in the calculations is given by 


W(k) = 


2ne 2 

k( 1 + 27TX2Dk) 


(33) 


with ^ 2 d = 6.6 A for intrinsic M 0 S 2 . Results are presented 
for a 120 x 120 sampling of the Brillouin zone, which we have 
found necessary to converge the binding energy to roughly 0.1 
eV accuracy, in agreement with the fully ab initio BSE study 
presented in Ref. [7]. Specifically, for M 0 S 2 this sampling 
gives a E exciton binding energy of 0.41 eV, however an ex¬ 
trapolation to the infinite sampling limit gives approximately 
0.52 eV, in good agreement with our prior result obtained in 
Ref. [6] (0.54 eV). In Fig. 4, the conduction bands have been 
rigidly shifted to increase the band gap to 2.41 eV, such that 
the Is exciton peak occurs near its experimentally observed 
value of 2.0 eV (due to the spin-orbit interaction, this peak is 
actually split into the so-called A and B peaks at about 1.9 and 
2.0 eV respectively [2]). An important conclusion to be drawn 
from Fig. 4 is that the more realistic band structure generates 
only minor quantitative differences in £ 2 ( 01 ), compared to that 
generated by the two band model. 

The labeling of states in Fig. 4 is done via inspection 
of the wavefunction, in either reciprocal or real-space. For 
example, in Fig. 5 we show the selection-rule-determining 
product A* c (k)P v _?(k) [which is closely related to the pseudo- 
wavefunction A* (&)] for right-handed polarization. The sym¬ 
metries of the exciton wavefunctions are apparent, and the val¬ 
ley selectivity is also recovered in the presence of excitonic 
effects. 

Focusing on the features in the £ 2 (oj) spectrum that derive 
from the s-type exciton states, the Rydberg series is nonhy- 
drogenic, as discussed in detail in Refs. [10, 16]. This follows 
from the unusual form of the screened Coulomb interaction 
for these monolayer thick materials. In particular, it deviates 
substantially from the 1 /s^r form that dominates in conven¬ 
tional semiconductors. The Hamiltonian with this latter inter¬ 
action has additional symmetry which leads to the “acciden¬ 
tal” angular momentum degeneracy in the hydrogen spectrum. 
Here that symmetry is broken: we find that for a given princi¬ 
pal quantum number, the larger angular momentum states are 
more strongly bound, i.e. E\ s < E 2 P < E 2 s < £ 3 d and so 
on. The same behavior has been recently observed in a fully 
ab initio BSE calculation [19], and the present work provides 
a simple physical explanation for this behavior in terms of 
the effective screened interaction (see also Refs. [33, 34] for 
similar findings). To verify this unconventional disposition 
of dark exciton states requires a nonlinear spectroscopic mea¬ 
surement, which we discuss in the next section. Furthermore, 
we also note a small splitting of the 2 p, 3d , and 3 p dark exci¬ 
ton states. In particular, the 20 meV splitting of the 2 p states is 
in good agreement with recent results [31, 33]. As mentioned 
before, Srivastava and Imamoglu have traced this degeneracy 
breaking to the orbital overlaps in Eq. (22) and explained the 
effect in terms of Berry curvature in the single-particle band 
structure [31]. 


IV. TWO-PHOTON ABSORPTION 

Our theoretical framework for the two-photon absorption 
essentially follows the early work of Mahan [35] for 3D semi¬ 
conductors and Shimizu [36] for 2D quantum wells including 
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Slater determinant and that p is diagonal in reciprocal space. 
To have a nonzero Eq. (35) requires that the real-space exci- 
ton wavefunctions A F and A M have orbital angular momenta 
which differ by ± 1; this is the same two-photon selection rule 
as found in conventional semiconductors including consider¬ 
ation of exciton effects. Combined with the result of the pre¬ 
vious section - that one-photon absorption produces s-type 
excitons - we conclude that two-photon absorption produces 
only p-type excitons. With these results, the two-photon ab¬ 
sorption essentially follows the early work of Mahan [35] for 
3D semiconductors or Shimizu [36] for 2D quantum wells. 

The primary complication in the evaluation of two-photon 
absorption is the evaluation of the internal sum over interme¬ 
diate states in Eq. (34). We follow the approximation intro¬ 
duced by Mahan [35] and used by Shimizu [36] that allows the 
sum to be eliminated with a completeness relation. Explicitly 
incorporating the above results, the first term in Eq. (34) can 
be written as (the second term is analogous) 



A%(r = 0)Aff»(r) _ 

Em — Ej — fuO] 


VrA£(r) 


-itiPo 

Eg - (E b ) - hcoi 


[V-v.aJ»L =o 


(36) 


FIG. 5. Reciprocal space plots of the selection-rule-determining 
product Af c (k)P™(k). In the presence of right-handed circular po¬ 
larization, it is seen that excitons are only created at the K point, and 
not at the K' point, as was found in Ref. [3] and in Sec. IIC in the 
absence of exciton effects. 


explicit consideration of excitons. For a two-photon process, 
the transition rate is again given by Eq. (11), except we now 
have two perturbing fields, V* = ( eAi/mc)Xi • p (i = 1,2), 
where A t is the vector potential, A; is the polarization, and ficoi 
is the photon energy. The matrix element of the perturbation 
can be evaluated by a sum over intermediate states |M), 


V/f= (£) 


</|A 2 • p|M)(M|Ai • p\F) 

Em — Ej — Ti(jj\ 
{I\\ x -p\M){M\\ 2 -p\F) 


Em — Ej — haj 2 


(34) 


The two-photon spectroscopy of single-particle states is triv¬ 
ial, and so we restrict our analysis to the excitonic case. As in 
the one-photon exciton absorption, Eq. (31) holds for the ma¬ 
trix element connecting the ground and intermediate exciton 
states. In contrast, the matrix element between two exciton 
states (intermediate and final) is 

<M|Aj -p\F) = fc^Af (fc)Ai • kA F vc (k) 

k 

= nj] Af (fc)e- ,T0fc Ai • kA F vc (k)e h ^ (35) 

k 

= -ih J d 2 rAf c *(r) A! • V r A F c (r). 

In the above, we have restricted the analysis to two bands (c, v) 
and used the facts that the expectation value of p is zero in a 


where (Eb) is an average intermediate (s-type) exciton energy 
introduced to facilitate the (complete) sum over intermediate 
states; for simplicity we will henceforth set (Eb) to zero as its 
primary influence is to simply alter the prefactor. In contrast 
to the hydrogenic exciton case, where further results can be 
obtained analytically, the matrix elements here must be eval¬ 
uated numerically. 

The two-photon transition rate is thus given by 


W(n) = 2nh (—— ) 4 (A!A 2 ) 2 (HPo) 2 
\mcJ 


z 


[A, • V,. i'!n| r 


Ea ~ fl(j)\ 


+ {1 ~ 2} 


5(h£l - E F ) 

(37) 


where TiEl = ti<j)i+ h(L> 2 . The simplest case to consider is when 
Ai = A 2 and ?ko\ = fico 2 ~ E g / 2, which gives 

MQ) = Wo E I A • VrA F c (r)\l = 0 6(hn - Ep) (38) 

F 


where 


Wo 


327zTz 3 e A A 2 x A\P^ 
m A c A E 2 g 


(39) 


If both photons have the same circular polarization, then this 
experiment probes valley-selective p- type excitons, which are 
dark in the linear measurement. Using photons with opposite 
polarizations would create /7-type excitons in both valleys. 

Motivated by recent nonlinear spectroscopic measurements 
on WSe 2 [9] and WS 2 [19], in Fig. 6 we show the results 
of a numerical evaluation of Eq. (38) for these two materi¬ 
als; the exciton wavefunctions and their derivatives have been 
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obtained from the real-space effective mass treatment of the 
two-band model (i.e. the small splitting of the p excitons is ne¬ 
glected). The agreement with experiment, for both the linear 
and nonlinear response, is seen to be quite good. In the calcu¬ 
lations, we have used the same screening length, xid - 7.0 A 
for both materials, which yields an exciton binding energy of 
0.48 eV (in accord with our previous results [6]). We note 
that this exciton binding energy is slightly larger than that de¬ 
termined in Refs. [9, 10] (0.37 and 0.32 eV for WSe 2 and WS 2 
respectively). 

In the narrow linewidth limit, the two-photon absorption 
identifies the /7-type excitons with energies slightly below that 
of the corresponding s-type exciton. For a larger linewidth, 
the 2 p transition is still resolved and responsible for the 
main peak seen in experiment, while the remaining transi¬ 
tions merge to yield a weak feature before the continuum on¬ 
set. Importantly, ratio between the 2 p peak height and the 
higher-energy signal (near the continuum onset) is determined 
by the spectral linewidth. It is thus encouraging that our sim¬ 
ulated spectrum simultaneously reproduces the 2 p linewidth 
and this intensity ratio; the required broadening suggests that 
it should be difficult to observe the 3 p transition at this resolu¬ 
tion. This leaves open the origin of the small feature observed 
near 2.5 eV in the experimental spectrum for WS 2 . 

Finally, we point out that a recent study on WSe 2 us¬ 
ing one- and two-photon photoluminescence excitation spec¬ 
troscopy [20], has identified the 2s and 2 p transitions to have 
the same energy to meV accuracy. This is in quite stark con¬ 
trast with the results of the present work, which suggest that 
the 2 p exciton should be lower in energy by at least 50 meV. 
We hope that future work, both experimental and theoretical, 
is devoted to investigating this discrepancy. 


V. CONCLUSIONS 

In this work, we have expanded the effective mass the¬ 
ory presented in Refs. [6, 10] to include a fully k-dependent 
model of the band structure, in harmony with other recent 
works [8, 33, 37]. This extension allows for deviations from 
parabolicity, including trigonal warping behavior which has 
been emphasized in other contexts [23, 38]. We find that two- 
and three-band models of the single-particle band structure 
give nearly identical results for the exciton properties within 
a simplified BSE formalism, suggesting that trigonal warping 
is a secondary effect. Furthermore, our numerical results are 
nearly identical to those of the effective mass treatment from 
our previous work [6, 10], justifying its use in those contexts. 
We have definitively proved that spin-singlet excitons with 
s-type azimuthal symmetry, which have been the most stud¬ 
ied [6, 8, 10], are indeed the optically bright excitons. As in 
our previous work [10, 16], we confirm that the disposition of 
bright exciton states is distinctly non-hydrogenic. 

The dark spin-singlet excitons have also been investigated 
and found to exhibit another deviation from the hydrogen 
model, in the form of a broken angular momentum degen¬ 
eracy. Using an approach similar to ours, the authors of 
Refs. [33, 34] have identified the same qualitative behavior. 


Abs theory - TP A theory 

• Abs expt • TPPLE expt 



FIG. 6. Two-photon absorption (TPA) intensity for monolayer 
(a) WSe 2 and (b) WS 2 evaluated numerically with Eq. (38) (blue 
lines). The spectra have been artificially broadened with a Gaus¬ 
sian linewidth (FWHM) of 80 meV (thicker line) and 20 meV (thin¬ 
ner line). The experimental two-photon photoluminescence excita¬ 
tion (TPPLE) spectrum for WSe 2 [9] and WS 2 [19] is included for 
comparison (blue circles). The theoretical linear absorption spec¬ 
trum from the same model (FWHM of 50 meV) is overlaid for refer¬ 
ence (grey lines) along with the experimental result (gray circles) for 
WSe 2 [9] andWS 2 [15]. 


This observation will be key in future analyses of two-photon 
spectroscopies on TMDCs. A recent manuscript contains re¬ 
sults from a fully ab initio BSE calculation on WS 2 and also 
finds this peculiar angular momentum behavior [19]. It is 
clearly encouraging that our simple low-energy theory - fea¬ 
turing a few-band representation of the single-particle states 
and an appropriate treatment of screening with a model di¬ 
electric function - is able to correctly reproduce the optical 
selection rules, the character of bright and dark exciton states, 
the broken angular momentum degeneracy, the quantitatively 
large exciton binding energies, and the spectral features of the 
nonlinear two-photon absorption. In this regard, we believe 
the model presented here represents perhaps the simplest pre¬ 
dictive minimal model capable of unifying these wide-ranging 
features in monolayer TMDCs. 

Note added. As discussed in the main text, a recent preprint 
analyzes the impact of the band overlap factors in the effective 
Coulomb interaction, Eq. (22), and systematically develops 
the next order terms in k-k', demonstrating signatures of the 
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Berry curvature in the exciton spectra [31]. Our numerical 
results agree with their analysis and with their estimate for the 
splitting of the 2 p exciton levels. Figure 4 was updated to 
reflect these splittings. 
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